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Abstract 

A C implementation of Kruskal’s up-and-down-blocks monotone regression algorithm 
for use wth .C(), and a comparison with other implementations. 
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Note: This is a working paper which will be expanded/updated frequently. All suggestions 
for improvement are welcome. The directory delceuwpdx.net/pubfolders/jbkPava has a pdf 
version, the bib hie, the complete Rmd hie with the code chunks, and the R source code. 

1 Introduction 

O no ! Not another monotone regression implementation !! There are already so many !!! 

There is isoregO in the stats package (R Development. Core Team (2017)), gpavaO in 
isotone() (De Leeuw, Hornik, and Mair (2009)), and pava in Iso (Turner (2015)). There 
is also wmonregO in smacof (De Leeuw and Mair (2009), Mail', De Leeuw, and Groenen 
(2015)) and amalgmO from De Leeuw (2016). Now gpavaO is written in R, amalgmO calls 
the Fortran from Cra.n (1980), pava() from Iso calls ratfor Fortran, and isoregO only 
does unweighted least squares. I wanted something in C (because C is not Fortran) which 
did weighted monotone regression, and I wanted to use the .C() interface (because I am 
exceedingly simple). 
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Now wmonreg from smacof fits the bill. It was written in 2014 by Patrick Groenen and 
Gertjan van den Burg. Same as our proposed algorithm here, it implements the up-and-down- 
blocks algorithm of Kruskal (1964). If we compare computation time then wmonreg() is the 
main competitor. But I also wanted code that was easy to modify for different unimodal loss 
functions and that performed relatively uniformly over the range of “almost in the correct 
order” to “order completely wrong”. 


2 Algorithm 

The best possible description of the algorithm was already given by Kruskal (1964) (page 
127). We copy his recipe, with a slight change of notation and terminology. 

“Our algorithm starts with the finest possible partitions into blocks, and joins the blocks 
together step by step until the correct partition is found. The finest possible partition consists 
naturally of n blocks, each containing only a single re;. 

Suppose we have any partition into consecutive blocks. We shall use Xb to denote the average 
of the Xi in block b. If 6_, b , b + are three adjacent blocks in ascending order, then we call b 
up-satisfied if Xb < %b + and down- satisfied if Xb _ < Xb ■ We also call b up-satisfied if it is the 
highest block, and down-satisfied if it is the lowest block. 

At each stage of the algorithm we have a partition into blocks. Furthermore, one of these 
blocks is active. The active block may be up-active or down-active. At the beginning, 
the lowest block, consisting of x m i n , is up-active. The algorithm proceeds as follows. If 
the active block is up-active, check to see whether it is up-satisfied. If it is, the partition 
remains unchanged but the active block becomes down-active; if not, the active block is joined 
with the next higher block, thus changing the partition, and the new larger block becomes 
down-active. On the other hand, if the active block is down-active, do the same thing but 
upside-down. In other words, check to see whether the down-active block is down-satisfied. 
If it is, the partition remains unchanged but the active block becomes up-active; if not, the 
active block is joined with the next lower block into a new block which becomes up-active. 
Eventually this alternation between up-active and down- active results in an active block 
which is simultaneously up-satisfied and down-satisfied. When this happens, no further 
joinings can occur by this procedure, and we transfer activity up to the next higher block, 
which becomes up-active. The alternation is again performed until a block results which is 
simultaneously up-satisfied and down-satisfied. Activity is then again transferred to the next 
higher block, and so forth until the highest block is up-satisfied and down-satisfied. Then the 
algorithm is finished and the correct partition has been obtained." 

Our implementation jbkPavaO, which uses Kruskal’s initials to distinguish it from other 
pava’s, stays as close as possible to this description. Blocks are C structures, collected in an 
array. Each block has a value, a weight, a size, the index of the previous block, and the index 
of the the next block. 

struct block { 
double value; 
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double weight; 
int size; 
int previous; 
int next; 

>; 

At the end the block values are expanded and returned in the original vector. Space for 
the blocks is allocated and freed in the routine (i.e. it is not under the control of R). Given 
Kruskal’s description of the up-and-down blocks algorithm, the code in the appendix should 
be pretty readable. 


3 Timings 

We apply the six monotone regression methods to vectors of length 10,000 with unit weights, 
repeating each analysis 100 times. We report the elapsed time from system.time(). 


3.1 Random 

The vectors of length 10,000 are normal random numbers (a different vector for each of the 


runs). 


## 

wmonreg 

0.386 

## 

gpava 

280.351 

## 

isoreg 

0.128 

## 

amalgm 

2.550 

## 

Iso::pava 

6.544 

## 

jbkPava 

0.137 


Clearly jbkPavaO and isoregO are the clear winners, although isoregO has the advantage 
that it does not have to take weighted means and keep track of the block weights. Note that 
gpavaO, written in R, is abysmally slow. The two methods amalgmO and Iso: :pava() 
that call Fortran routines are not really competitive. wmonregO is doing quite well, but 
jbkPavaO is more than twice as fast. 


3.2 Reversed 


We apply our six methods, one hundred times each, to the vector 10000:1, which obviously 
will produce a completely merged vector with all elements equal to the mean. 


## wmonreg 
## gpava 
## isoreg 


4.424 

425.842 


0.027 
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## amalgm 2.373 

## Iso::pava 10.727 

## jbkPava 0.052 

In this situtation, in which blocks are never up-satisfied and always down-satisfied, again 
jbkPavaO and isoregO are best. wmonregO and Iso: :pava() lose ground. 

3.3 Ordered 

We apply our six methods, one hundred times each, to the vector 1:10000, which obviously 
just returns the vector itself. There is no merging at all, block are always up-satisfied and 
down-satisfied, and we merely loop through the vector. 


## 

wmonreg 

0.026 

## 

gpava 

5.896 

## 

isoreg 

20.162 

## 

amalgm 

0.049 

## 

Iso::pava 

0.011 

## 

jbkPava 

0.042 


wmonregO is best-in-show, which is interesting because in many applications we expect the 
vector to be closer to the correct ordering than to random or reverse ordering. Both Fortran 
calling routines perform well. isoregO falls flat on its face and is even beaten by gpava(). 
There is probably some reasonable explanation for this, but since unweighted isoregO is 
not really in the competition I did not explore this. 

3.4 Conclusion 

Both pavaO and wmonregO perform well over the range of examples. There is an indication 
that pavaO is considerably better for vectors that are out of order, while wmonregO may be 
better for vectors that are already close to correct. In future versions of this paper we will 
try to speed up pavaO performance. 

4 Code 

4.1 R code 

dyn.load(" jbkPava. so") 

jbkPava <- function (x, w = rep(l, length(x))) { 
h <- 

.C(" jbkPava" , 
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x = as.double (x), 
w = as.double (w), 
n = as.integer(length (x))) 
return (h$x) 

> 

4.2 C code 

#include <stdlib.h> 

#include <stdbool.h> 

struct block { 
double value; 
double weight; 
int size; 
int previous; 
int next; 

>; 

void jbkPava (double *x, double *w, const int *n) { 

struct block ^blocks = calloc ((size_t) * n, sizeof(struct block)); 
for (int i = 0; i < *n; i++) { 
blocks[i].value = x[i]; 
blocks[i].weight = w[i]; 
blocks[i].size = 1; 
blocks [i].previous = i - 1; 
blocks[i].next = i + 1; 

> 

int active = 0; 
do { 

bool upsatisfied = false; 
int next = blocks[active].next; 
if (next == *n) upsatisfied = true; 
else if (blocks[next].value > blocks[active].value) 
upsatisfied = true; 
if (!upsatisfied) { 

double ww = blocks [active].weight + blocks[next].weight; 
int nextnext = blocks[next].next; 

double wxactive = blocks[active].weight * blocks[active].value; 

double wxnext = blocks[next].weight * blocks[next].value; 

blocks[active].value = (wxactive + wxnext) / ww; 

blocks[active].weight = ww; 

blocks[active].size += blocks[next].size; 

blocks[active].next = nextnext; 
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if (nextnext < *n) 

blocks[nextnext].previous = active; 
blocks[next].size = 0; 

} 

bool downsatisfied = false; 
int previous = blocks[active].previous; 
if (previous == -1) downsatisfied = true; 
else if (blocks[previous].value < blocks[active].value) 
downsatisfied = true; 
if (!downsatisfied) { 

double ww = blocks [active].weight + blocks [previous].weight; 

int previousprevious = blocks[previous].previous; 

double wxactive = blocks[active].weight * blocks[active].value; 

double wxprevious = blocks [previous].weight * blocks[previous].value; 

blocks[active].value = (wxactive + wxprevious) / ww; 

blocks[active].weight = ww; 

blocks[active].size += blocks[previous].size; 
blocks[active].previous = previousprevious; 
if (previousprevious > -1) 

blocks[previousprevious].next = active; 
blocks[previous].size = 0; 

} 

if ((blocks [active].next == *n) && downsatisfied) break; 
if (upsatisfied && downsatisfied) active = next; 

} while (true); 
int k = 0; 

for (int i = 0; i < *n; i++) { 
int blksize = blocks[i].size; 
if (blksize >0.0) { 

for (int j = 0; j < blksize; j++) { 
x[k] = blocks[i].value; 
k++; 

} 

} 

> 

free (blocks); 

> 
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